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We present a new Monte Carlo muon propagation algorithm MUM (MUons+Medium) which 
possesses some advantages over analogous algorithms presently in use. The most important features 
of algorithm are described. Results on the test for accuracy of treatment the muon energy loss 
with MUM are presented and analyzed. It is evaluated to be of 2xl0 -3 or better, depending 
upon simulation parameters. Contributions of different simplifications which are applied at Monte 
Carlo muon transportation to the resulting error are considered and ranked. It is shown that when 
simulating muon propagation through medium it is quite enough to account only for fluctuations 
in radiative energy loss with fraction of energy lost being as large as 0.05-^0.1. Selected results 
obtained with MUM are given and compared with ones from other algorithms. 
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I. INTRODUCTION 



Muon propagation through thick layers of matter was in the scope of interest for a long time since the first under- 
ground experiments with natural muon and neutrino fluxes had started. Development of "underground" technique has 
led to creation of number of underground, underwater and under-ice detectors by which a wide spectrum of problems 
is presently under investigation. Accurate calculation of muon transport plays an important role for such experiments 
because (a) neutrinos are detected by muons which are born in vN interactions and propagate a distance in medium 
from the point of interaction to a detector; (b) muons which are produced in atmospheric showers generated by cosmic 
rays represent the principal background for neutrino signal and therefore their flux at large depths should be well 
known; (c) atmospheric muons deep under sea or earth surface are the only intensive and more or less known natural 
calibration source which allows to confirm correctness of the detector model by comparison experimental and expected 
detector response; (d) flux of atmospheric muons itself carries the physical information which is of interest. 

Along with analytical and semi-analytical methods (Refs. [p]-p^|) one widely uses Monte Carlo (MC) technique 
(Refs. [|14H25[) which directly accounts for the stochastic nature of muon energy losses to simulate the muon propa- 
gation through matter. There are several MC muon transportation algorithms currently in use (see, e.g. Ref. ]26[ ] for 
detailed analysis of their advantages and disadvantages), but essential theoretical and experimental progress of last 
years makes to create new ones. Here we present a MC muon propagation code MUM (MUons+Medium) written 
in FORTRAN which possesses some advantages in accuracy and flexibility over analogous simulation algorithms (al- 
though it does not contain some important features in its current version, e.g. it does not give the 3D information 
about angular and lateral deviations of muons). The algorithm has been developed for the Baikal deep underwater 
neutrino experiment (Ref. |27j ) but we believe it to be useful also for other experiments with natural fluxes of high 
energy muons and neutrinos. When working on MUM we aimed at creation of an algorithm which would: 

(a) account for the most recent corrections for the muon cross sections; 

(b) be of adequate and known accuracy, i.e. does not contribute an additional systematic error which would exceed 
one from "insurmountable" uncertainties (e.g. with muon and neutrino spectra and cross sections) and whose 
value would be well known for any setting of simulation parameters; 

(c) be flexible enough, i.e. could be optimized for each concrete purpose to desirable and well understood equilibrium 
between computation time and accuracy and be easily extended for any medium and any correction for the cross 
sections of the processes in which high energy muon looses its energy; 

(d) be "transparent" , i.e. provide user with the whole set of data related to used models for the muon cross sections; 

(e) be as fast as possible. 

We describe the main features of our algorithm in Sec. H. Sec. [II gives an analysis for the algorithm accuracy. In 
Sec. IV we report the results of investigation on the set of parameters which should be used to simulate the muon 



1 



propagation with an optimum equilibrium between accuracy and computation time. Sec. ^ presents selected results 
obtained with MUM in comparison with ones from other muon propagation MC codes, namely PROPMU (Ref. [ p2| ) 
and MUSIC (Ref. p3). Sec. VI gives general conclusions. We also present parameterizations for muon cross sections 



as they are used in MUM in Appendix J 
treated in our algorithm, in Appendix [E 



and give proof of formula for free path between two muon interactions, as 



II. ALGORITHM DESCRIPTION 



The basic features of the MUM algorithm are as follows. 

(a) The code does not use any preliminary computed data as an input, all necessary tables are prepared at the 
stage of initiation on the base of five relatively short routines, four of which return differential cross sections 
da(E,v)/dv (where E is muon energy, v is fraction of energy lost v = AE/E) for bremsstrahlung, direct e + e~- 
pair production, photo-nuclear interaction and knock-on electron production, correspondingly, and fifth one 
does stopping power due to ionization [dE(E)/dx]i on (see Appendix [a| for corresponding formulas). Thus, it is 
easy for any user to correct or even entirely change the model for muon interactions as it is necessary. Also any 
material can be easily composed. Three media, namely pure water, ice and standard rock are available directly. 

(b) We have tried to decrease the "methodical" part of systematic error which originates from finite accuracy of 
numerical procedures on interpolation, integration, etc., down to as low level as possible, special attention was 
put on procedures which simulate free path between two sequential muon interactions and fraction of energy 
lost. To combine this with the high speed of simulation the values for free paths, energy losses, differential and 
total cross sections along with solutions for all ordinary and integral equations are computed in MUM at the 
initiation stage, tabulated and then referenced when necessary with an interpolation algorithm whose accuracy 
has been carefully tested for each table by comparison with directly computed values to be not worse than 0.5% 
(typically, much better). 

(c) The most important parameters are changeable and can be tuned to an optimum combination, depending upon 
desirable accuracy, necessary statistics and restriction on the computation time for each concrete problem. 

(d) The code combines algorithms for muon transportation through thick layers of matter down to detector and 
for simulation muon interactions within detector sensitive volume (these algorithms have to differ from each 
other). This is important for deep underwater and under- ice Cherenkov neutrino telescopes (see Refs. [^7| j30||), 
where the same material (water or ice) represents both a shield which absorbs atmospheric muons and detecting 
medium in which muons and shower particles resulting from muon interactions generate Cherenkov photons 
detected by phototubes. 

(e) Formally, initial muon energies up to 1 EeV can be processed by the MUM algorithm but uncertainties with 
muon cross sections which grows along with the muon energy (especially, for photo-nuclear interaction) make 
one to apply MUM output with care at muon energies E > Q.l-j-1 PeV. Landau-Pomeranchuk-Migdal effect is 
not accounted in MUM. 

(d) Besides muon transportation algorithm itself, the code includes number of routines which allows to obtain 
directly values for differential and total cross sections, mean free paths, energy losses and other related data 
for the given set of input parameters. Sampling the atmospheric muon energies at the sea level according to 
different models for spectrum is possible with MUM, as well. Also several test procedures are included which 
provide by data concerning accuracy of different algorithm steps. See Sec. Ill, Sec. |y| and Sec. |v| for selected 
output of these procedures. 

The usual approximation for treatment the muon energy loss is applied in the MUM algorithm: muon interactions 
with comparatively large energy transfers when fraction of energy lost v exceeds some value v cu t, are accounted 
by direct simulation of v > v cut for each single interaction according to shape of differential cross sections (these 
interactions lead to "stochastic" energy loss or SEL) while the part of interaction with relatively small v is treated by 
the approximate concept of "continuous" energy loss (CEL) using the stopping power formula 
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Here index j indicates type of interaction (J = b for bremsstrahlung, j = p for direct e + e~-pair production, j = n 
for photo- nuclear interaction and j = e for knock-on electron production, respectively); index i runs over n kinds of 
atoms given material consists of; fcj = Ni/N to t is fraction of i-th element; N% and N to t are number of given kind of 
atoms and total number of atoms, respectively, per unit of material volume; Na is the Avogadro number; p is the 
material density; A e tf = N^ ot Y^i=i(NiAi) is an effective atomic weight for given material; A4 is atomic weight for 
i-th element; v^ in is minimum kinematically allowed fraction of energy lost for i-th element at j-th process. One is 
forced to decompose energy losses into two parts because simulation of all interactions with v > v m i n would result 
in infinite computation time due to steep dependence of muon cross sections on v: they decrease with v at least as 
da(E,v)/dv oc v and for some processes are not finite at v — > 0. Number of interactions to be simulated per unit 
of muon path grows, roughly, as N int oc v~ ut along with computation time. Actually, two different criteria by which 
the given muon interaction is attributed either to SEL or to CEL are available in the frame of the MUM algorithm. 
The first one (relative) has been described above and is applied when muon is transported down to detector location. 
Second, absolute criterium, is useful when simulating muon interactions within an underwater or under-ice array to 
obtain the detector response with the fixed energy threshold: interaction is of SEL type if AE > AE cu t and of CEL 
type, if AE < AE cut . Optionally, cross section for knock-on electron production da e (E,v)/dv can be set in MUM 
to zero, in which case muon propagation down to detector is simulated with entirely "continuous" ionization and 
its fluctuations are neglected but simulation of the muon interactions with fixed energy threshold includes knock-on 
electron production in any case. Both v cut and AE cut represent parameters for MUM initiation procedure and can be 
set to any values within 1CU 4 < v cut < 0.2 and 10 MeV < AE cut < 500 MeV, correspondingly. The optimum value for 
AE cut depends upon configuration of the given detector and also upon characteristics of algorithms which simulate 
the shower development, the Cherenkov photons generation and propagation, and detector response. All this is ou t 
of gi ven a rticle scope, therefore we discuss this parameter nowhere below except for mentioning that Eqs. (2T), (2.3) 
and (2.6) are used in algorithm with absolute treatment of the muon energy loss decomposition being modified by 
replacement v cu t — > AE cu t/E. Influence of simplified entirely "continuous" treatment of ionization and value of v c ut 
upon simulation accuracy is analyzed in details in Sec. [V (see also Ref. J3l|). 
The principal steps of simulation are as follows. 

/. For muon with initial energy E\ the free path L till interaction with v > v cut is simulated. For this, after a 
random number 77 uniformly distributed in a range from to 1 has been sampled, one solves the following set of 
equations: 
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(proof of Eqs. ( |2.2|) is given in Appendix |B|). Here, E% < E\ is muon energy at the point of interaction and energy 
dependent mean free path L(E) between two interactions with fraction of energy lost v > v cut is expressed by 
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where v\^ ax is maximum kinematically allowed fraction of energy lost for i-th element at j'-th kind of interaction. 
First equation in Eqs. (2.2) is solved for variable E2, then free path L can be found by second equation. We would 



like to stress that such approach allows to perform the accurate simulation independently on chosen value of v cu t in 
contrast to commonly used simplification 
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which neglects dependence L(E) upon energy and, consequently is (a) the less accurate the larger v cu t is and (b) 
produces the error of different signs for the cases when ionization is included in SEL or its fluctuations are neglected. 
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It is illustrated by two plots in Fig. [j]. Upper plot shows function L(E) for pure water. Two sets of curves are 
presented for two models of ionization. Each set includes dependencies for 3 values of v cu t. 1CP 4 , 1CP 3 and 1CP 2 . 
In fact, L(E) is almost a constant at E > 5 TeV but changes steeply at lower energies. It increases with decrease 
of energy if ionization is entirely "continuous" and, on the contrary, it decreases if ionization is included in SEL. 



Thus, simulating free path by Eq. (2.4) one overestimates it (and consequently underestimates energy loss) in case 
when ionization is included in SEL and, on the contrary one underestimates free path and overestimates energy loss 
in case of completely "continuous" ionization. Lower plot in Fig. |l| shows resulting error in the value of simulated 
free path if — ln(r/) = 1 (for larger — ln(7y) the effect is more significant). The set of curves represent dependencies 



k(E) = L approx {E) j 'L(E) with L(E) computed by Eqs. ( [2.2; ) and L approx (E) computed by Eq. (2.4). With ionization 
included in SEL overestimation for free path is less than 1% at v cut < but reaches ~15% at v cut = 1CP 2 which 
leads to l-j-2% underestimation of total energy loss below muon energy 1 TeV. In case with "continuous" ionization 
the effect is of opposite sign and again is more significant for large v cut . 

II. After free path L and muon energy E2 have been found from Eqs. (2.2), the type of interaction is simulated 
according to proportion between total cross sections of different processes: 
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which are computed as: 
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III. Fraction of energy lost v is simulated according to shape of differential cross section for given process j: 
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(2.7) 



and new muon energy E 1 = E% ■ (1 — v) is determined. 

IV. Steps 7-777 are repeated sequentially until muon either reaches the level of observation or stops. Muon are 
considered as stopped as soon as its energy decreases down to 0.16 GeV which corresponds to the Cherenkov threshold 
for muon in pure water. 
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(a) - mean free path between two sequential muon interactions with fraction of energy lost v > v cu t 



FIG. 1. (a) - mean free path between two sequential muon interactions with fraction of energy lost v > v cu t L(E) (Eq. (2.3)) 
in pure water vs. muon energy. Two sets of curves correspond to two models of ionization. Thick lines are for ionization included 
in SEL, thin ones correspond to entirely "continuous" ionization. Solid lines: v cut = 10~ 4 ; dashed lines: v cut = 10~ 3 ; dotted 
lines: v cu t = 10 -2 . (b) - Function k(E) — L avvrox (E) j L(E) for — 111(77) ~ 1 ( see text). Thickness and type of lines are of the 
same meaning as for plot (a). 
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III. ALGORITHM ACCURACY 



As described in Sec. || the MUM code (as well as any muon MC propagation algorithm) consists of the set of 
procedures on numerical solution of equations, interpolation and integration. All these procedures are of finite accuracy 
and, consequently, the incoming model for muon energy loss is somewhat corrupted by them. Thus, resulting energy 
loss as it simulated by a code are not the same as energy loss as it can be calculated by integration the differential cross 
sections which are at the input of the same code. The difference between simulated and calculated energy loss contains 
errors which are contributed by each step of simulation algorithm and thus, is a good quantitative criterium for its 
inner accuracy, whose contribution to the resulting error must not exceed one which comes, e.g. from uncertainties 
with muons cross sections and medium composition. Therefore to demonstrate accuracy of presented algorithm we 
have chosen just data on the relative difference (L s — Li) /Li between simulated L s and integrated Li total muon 
energy loss as was obtained with MUM for the pure water (Fig. |^) and standard rock (p = 2.65 g cm -3 , A = 22, Z 
= 11, Fig. ||). Inner accuracy is presented in the figure as a function of muon energy for several values of v cu t and 
two models of ionization energy loss. Values of L s were obtained as follows. For each muon energy E\ a distance D 
was chosen and propagation of N = 4xl0 6 muons over this distance was simulated. The condition N ■ D ^> L(E\) 
must be obeyed to obtain statistically significant result but, at the same time, D should be short enough to be passed 
by muons without decrease their energy down to zero, which practically leads to D = 0.5 m-^300 m depending upon 
muon energy, value of v cut and kind of medium. For each i-th muon its final energy E\ was fixed and then L s was 
calculated as 

L > = h(* (3 - 1} 

Li was computed as 

L i = hE 1 -E 2 ), (3.2) 
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FIG. 2. Relative difference (L s — Li)/Li between "simulated" L s (Eq. (3.1)) and "integrated" Li (Eqs. (3.2), (3.3)) total 
muon energy loss in the pure water. Horizontal solid line on each plot shows averaged over 24 tested muon energies value for 
(L s — Li) I Li which, additionally, is given at the upper left corner by the figure. Statistical error at lcr-level is shown at each 
point, (a) - ionization is included in SEL; (b) - ionization is entirely "continuous". 
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FIG. 3. The same as in Fig. | for standard rock (p = 2.65 g cm~ 3 , A = 22, Z = 11). 
where E2 was found as a solution of integral equation for the muon range: 
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Horizontal solid line on each plot shows the value for (L s — Li)/ Li averaged over 24 tested muon energies which, 
additionally, is given at upper left corner by a figure. Fig. |2|(a) and Fig. ||(a) indicate an excellent inner accuracy 
of the MUM algorithm with ionization included in SEL both for water and standard rock. Up to v cu t — 5xlCU 2 all 
points are within 0.6%-deviations which, besides, are of both signs, so averaged accuracy remains better than 10 -3 . 
Fig. ||(b) and Fig. ||(b) (which correspond to simplified completely "continuous" ionization) shows somewhat worse 
accuracy of the algorithm which falls down when v cut increases. Accuracy was found to be within 1% (with except 
for few points around muon energy E = 100 GeV) up to v cu t = 10~ 2 with averaged accuracy within 2xl0 -3 . This 
last value may be used as a conservative evaluation of inner accuracy for the MUM algorithm. Statistically significant 
likeness of plots obtained for water and standard rock can be seen. 

Thus, we conclude that assuming an optimistic evaluation of 1% for uncertainties in muon cross sections 
(Rcfs. 26 3^]) the inner inaccuracy of MUM does not exceed them for any v cu t < 5xl0 -2 if ionization is included in 
SEL and for any v cu t < 10 -2 if ionization is treated as entirely "continuous process" , independently upon material. 



IV. THE OPTIMUM SETTING OF SIMULATION PARAMETERS 

As was described in Sec. |0] v cu t is a parameter in the MUM algorithm and can be set optionally to different values. 
The larger is v cu t the higher is the speed of simulation, because the less muon interactions have to be simulated 
per unit of the muon path. But, on the other hand, too large value of v cut leads to the lost of accuracy since some 
essential part of fluctuations in the muon energy losses comes out of direct simulation. Thus, the question is how large 
value of v cu t may be chosen to keep result within desirable accuracy? Also different models for ionization can be used: 
it can be optionally either treated as completely "continuous" process or included in SEL. Small energy transfers 
strongly dominate at knock-on electron production (da e (E , v) / dv oc v~ 2 ), so this process is almost non-stochastic and 
it seems to be reasonable to exclude knock-on electrons from simulation procedure when simulating SEL which saves 
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computation time noticeable. How much does it affect the result o f si mulation? Influence of these factors on simulated 
result had been discussed in literature (see, e.g. Refs. |I^,^,^J,g5| ) but in our opinion more detailed analysis was 
lacking. Therefore we have undertaken our own investigation which is reported in this Section. For that we performed 
several sets of simulations both for propagation of mono-energetic muon beams and atmospheric muons sampled by 
sea level spectrum (in the later case we limited ourselves by simulation only vertical muons) through pure water down 
to depths from D = 1 km to D = 40 km. Of course, distances of more than several kilometers for vertical muons do 
not concern any real detector but simulations for large depths allow us to study general regularities which correspond, 
e.g. to nearly horizontal directions. Several runs were done for standard rock, as well. We tested different settings of 
parameters which were as follows. 

(a) v cu t, which changed within a range of 10 -4 < v cut < 0.2. Inner accuracy of the MUM code becomes somewhat 
worse at v cut > 5xl0 -2 , especially if fluctuations in ionization are not simulated (Sec. Ill) therefore results for 
v C ut — 0.1 and v cut — 0.2 are presented here only to illustrate some general qualitative regularities. 

(b) Model for ionization. 

(c) Parameterization for vertical sea level atmospheric muon spectrum. Two spectra were tested, namely one 
proposed in Ref. [B3| (basic): 
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and the Gaisser spectrum (Ref. [p4|): 
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(d) Parameterization for total cross section for absorption of a real photon by a nucleon at photo-nuclear interaction 



(t 7 jv which was treated both according to the Bezrukov-Bu gaev parameterization proposed in Ref. 35 (basic) 



and the ZEUS parameterization (Ref. |3q]) (see. Appendix A 2 for formulas) 



(e) A factor k a which all muon cross sections along with stopping power due to ionization were multiplied by to test 
influence of uncertainties in muon cross sections (and, consequently, in energy losses) upon result. We applied 
k a = 1.0 as a basic value but set also k a — 0.99 and k a = 1.01, which corresponds to decrease and increase of 
total energy loss by 1%, respectively. Note that it is an "optimistic" evaluation, the real accuracy of existing 
parameterization for muon cross sections is worse (see Refs. |26| , |3"2] ] ) . 

For each run we fixed the muon spectra at final and several interim depths. The differences between obtained 
spectra were a point of investigation. 

At the first set of simulations we propagated mono-energetic muon beams of 4 fixed initial energies E s = 1 TeV, 
10 TeV, 100 TeV and 10 PeV down to slant depths D = 3.2 km, 12 km, 23 km and 40 km, respectively, through 
pure water. The value D for each initial muon energy was chosen so that majority of muons had been stopped after 
propagation the given distance. This allows to track differences in simulated results obtained with different settings 
of parameters for all segments of muon beam path. In each case propagation of 10 6 muons was simulated. Fig. ^ 
shows resulting survival probabilities p = Nr)/N s (where N s = 10 6 is initial number of muons and Njj is number of 
muons which have survived after propagation down to the slant depth D) vs. v cu t for final and five interim values of 
D. Two curves are given on each plot for two models of ionization. Also results for k a — 1.00 ± 0.01 and for cr 7 jv 
parameterization according to Ref. ( f3q| ) are presented as simulated with the most accurate value v cu t — 10~ 4 . 

The following conclusions can be done. 

(a) In most cases except for some plots of the lower row and the left column in Fig. ^ (which corresponds to low 
survival probabilities and low muon initial energies, respectively) uncertainty in our knowledge of muon cross 
sections gives the principal effect which essentially exceeds ones from other tested parameters. 

(b) The difference between survival probabilities for two models of ionization is the less appreciable the larger muon 
energy is. It is quite understandable because at muon energies E < 1 TeV ionization represents the great bulk 
of total energy loss, and vice versa, it becomes minor at E > 1 TeV. Thus, contribution which is given by 
ionization at higher energies is small and, the more, its fluctuations do not play an important role. For muons 
with initial energies E 3> 1 TeV fluctuations in ionization become important only at very last part of muon 
path and "are not in time" to produce some noticeable effect. 
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(c) Generally, parameterizations for a 7 N as proposed in Refs. |3|,|3(| do not show a noticeable difference in terms 
of survival probabilities, in most cases it is within statistical error or exceeds it only slightly. 

(d) Increase of v cu t gives effect of both signs in survival probabilities: function p(v cu t) grows at the beginning of 
muon path and falls at the last part. The same "both-sign" dependencies are observed for ionization model. 

(e) For v cut < 0.05 there is almost no dependence of survival probability on v cut except for very last part of muon 
path where survival probability becomes small. Generally, dependence p(v cu t) is the less strong the larger initial 
muon energy is. 
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FIG. 4. Survival probabilities p = Nn/N s (where N s — 10 6 is initial number of muons in the beam and No is number of 
muons which have survived after propagation down to slant depth D in pure water) vs. v cut . Values of p were obtained as a 
result of simulation with MUM for mono-energetic muon beams with initial energies E s = 1 Tev (1st column of plots), 10 FeV 
(2nd column), 100 FeV (3rd column) and 10 PeV (4th column). Each column contains six plots which correspond to six slant 
depths D (which differs for different E 3 ). Closed circles represent survival probabilities which were simulated with ionization 
energy losses included in SEL, open ones correspond to computation with completely "continuous" model of ionization. Fwo 
horizontal solid lines on each plot show the value for survival probability computed with all muon cross sections multiplied by 
a factor k a = 1.01 (lower line) and fc CT = 0.99 (upper line) for v cu t = 10 -4 . Horizontal dotted lines correspond to v cut = 10~ 4 
and cross section for absorption of a real photon at photo-nuclear interaction parameterized according to Ref. |3(| instead of 
Ref. |}5| which is basic in MUM. Note different scales at Y-axis. 

The last item is illustrated complementary by Fig. ^ and Fig. || which show that for all initial energies E s simulated 
survival probability does not depend, in fact, on v cu t until 90% (for E s = 1 TeV) to 99.5% (for E s = 10 PeV) muons 
have been stopped. 
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D (km) D (km) 

FIG. 5. Survival probability p vs. slant depth D in pure water down to which propagation of mono-energetic muon beam 
with initial energy E s = 1 TeV (a), 10 TeV (b), 100 TeV (c) and 10 PeV (d) is simulated. On each plot 10 lettered curves 
which correspond to different values of v cu t are shown. Meaning of letters is as follows: A - v cu t = 10~ 4 , B - v cut = 2xl0~ 4 , 
C - v cut = 5xl0~ 4 , D - v cut = 10~ 3 , E - Vcut = 2xl0~ 3 , F - v cu t = 5xl0~ 3 , G - v cu t = 10~ 2 , H - v cu t = 2xl0~ 2 , I - v cut 
= 5xl0 -2 , J - Vcut = 10 _1 . This figure displays results which were obtained by simulation with ionization losses included in 
SEL. Statistical errors (which cause some un-smoothness of curves at small p) are not shown. 




FIG. 6. Relation p/po vs. p. p is survival probability for muon s of initial energy E s — 1 TeV (a), 10 TeV (b), 100 TeV (c) 
and 10 PeV (d) at propagation through pure water with ionization included in SEL as simulated for different values of v cu t', 
po is survival probability simulated under the same conditions for v cu t = 10~ 4 . Difference in p/po becomes noticeable only at 
p <10 _1 , i.e. at the very last part of muon beam path. 
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FIG. 7. Muon spectra resulting from 10 6 muons with initial energy E a — 1 TeV as simulated with four models for different 
slant depths D in pure water. The first three columns represent spectra obtained with ionization included in SEL for v cut = 
1CP 4 (1st column), 10 -2 (2nd column), and 0.2 (3rd column). 4th column contains spectra obtained for entirely "continuous" 
ionization and v cu t = 10 -4 . On each plot value of survival probability p is indicated without statistical error which does not 
exceed 1%. 



It was shown above what is result of simulations with different models of ionization and values of v cu t- It was a 
special point of interest for us to track how and why does it influence upon behavior of survival probability. Fig. |?j 
shows how muon spectrum resulting from mono-energetic muon beam with initial energy E s = 1 TeV transforms when 
its propagation being simulated through pure water down to the slant depth of 3.2 km. Results for four settings of 
parameters arc presented by four columns of plots. The first three columns represent spectra obtained with ionization 
included in SEL for v cu t = 10~ 4 , 1CP 2 and 0.2. 4th column contains spectra simulated with entirely "continuous" 
ionization and v cut = 10 -4 . The spectra grouped into the first column represent the most accurate tuning both for 
v cut and ionization model. The first three columns demonstrate that compactness of spectra at the same slant depth 
is the higher the more value of v cut is. Put your attention to the right edge of spectra which shifts toward low energies 
when v cu t increases (it is the most noticeably for v cu t = 0.2). The reason is that at any slant depth energy of the 
most energetic muons in simulated beam is determined by CEL. These muons due to statistical fluctuations did not 
undergo interactions with v > v cu t and, consequently, lost energy only by CEL which grows when v cu t increases. 
That is why the maximum energy in simulated muon beam is lower for large values of v cu t ■ Fraction of muons which 
did not undergo an "catastrophic" act with v > v cu t till given slant depth grows with increase of v cut because free 
path between two sequential interactions with v > v cu t grows approximately as L oc v cu t- It leads, in particular, to 
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distinctly visible separated picks in spectra for v cu t = 0.2 consisted just of muons which lost energy only by CEL. 
Also some deficit of low energy muons appears if one sets v cu t to a large value. In this case left edge of spectrum 
is provided only with muons which interacted with large fraction of energy lost while for smaller v cu t an additional 
fraction of muons comes here. As a result simulated spectrum of initial mono-energetic muons at given slant depth is 
more narrow if v cu t is large and, on the contrary, more wide if v cut is small. 

Now it is easy to understand how value of v cu t influences on simulated survival probabilities. When simulated muon 
beam goes through a medium loosing energy both in CEL and SEL processes, its spectrum is constantly shifting to 
the left (energy decreases). For v cu t — 10~ 4 the left part of spectrum reaches E = at a smaller slant depth comparing 
with larger v cu t and survival probability starts to decrease. At the same slant depth survival probability for v cu t = 
1CP 2 and v cut = 0.2 is still equal to 1. Thus, for the first part of path the survival probability is always larger for large 
v cut . At some slant depth (which is equal to ~2.8 km in the given case) compactness of spectra simulated with large 
v cu t starts to play an opposite role. Due to more powerful CEL muons stop faster comparing with accurate simulation. 
So at the final part of the beam path simulated survival probability for large v cu t decreases faster comparing with 
accurate simulation and, for instance, for v cut — 0.2 the rest of muon beam which reaches the slant depth D = 2.72 
km (37% of initial number of muons) completely vanishes within the next 30 m of path, while some fraction of muons 
simulated with v cut — 10~ 4 (0.07%) escapes down to the slant depth of D = 3.2 km. Qualitatively the same effect leads 
to the same consequences if one treats ionization as completely "continuous" process. Again, spectra becomes more 
narrow since fluctuations in ionization do not work and, as a consequence, survival probability becomes significantly 
higher comparing with simulation with accurate treatment of ionization at the beginning of muon beam path and 
falls down essentially faster at the final part of path. 

Results presented above showed the significant influence which both model of ionization and value of v cut have over 
survival probability for mono-energetic muon beam. But for practical purposes the more important is how these factors 
do work for real atmospheric muons with a power spectrum? In Fig. || we present intensity of vertical atmospheric 
muon flux / at different depths of pure water D from 1 km to 20 km vs. v cu t as simulated with muons sampled 



according to sea level spectrum Eq. (4.1). Simulation continued until 10 4 muons reached given depth. Curves for two 



models of ionization are shown for each depth along with results for k a — 1.00 ± 0.01 at v cu t = 10 , parameterization 



for (t 7 at from Ref. p6fl at v cu t = 10 , sea level muon spectrum Eq. (4.2) at v cu t = 10 and all energy losses treated 



entirely as CEL (for depths D < 5 km only). General conclusions for case with atmospheric muons are qualitatively 
the same as observed for mono-energetic muon beams but quantitatively the influence of v cu t and model of ionization 
energy losses on the resulting muon flux at large depths is much weaker. One can conclude the following: 

(a) Except for case D = 1 km computed muon flux is strongly affected by accounting for fluctuations in energy 
losses: muon flux intensity simulated with non-stochastic model of energy loss is less comparing with stochastic 
model by 10 % at 3 km w.e. and by 20% at 5 km w.e.. At the depth of 20 km of pure water vertical muon flux 
computed with ignorance of fluctuations is only 10 % of simulated flux. 

(b) Like a case with mono-energetic beams 1%-uncertainty in muon cross sections plays the principal role for resulting 
error in simulated muon depth intensity. This error has a tendency to grow with depth from ±2.5% at depth 
of 1 km w.e. to ~ ±15% at 20 km w.e.. But a particular case of this uncertainty, namely difference between 



parameterizations for cr 7 jv from Refs. 35 36|], does not lead to a significant difference in resulting intensity. 



(c) Difference between muon spectra Eq. ( |4.l| ) and Eq. fl4.2p leads to uncertainty from -4% (D = 1 km) to 16% (D 
= 20 km). 

(d) Error which appears due to simplified, entirely "continuous" ionization lies, commonly, at the level of 2-^3%. 

(e) Dependence of simulated muon flux intensity upon v cu t is the most weak one comparing with other studied error 
sources. Function I(v cu t) is almost a constant at v cu t < 0.05 and changes in a range ±1-^2% which is very close 
to statistical error. Up to v cut =0.1 the contributed error is less than one which comes from ± 1%-uncertainty 
with the muon cross sections. Also no statistically significant influence of v cu t upon the shape of differential 
atmospheric muon spectra was observed at all tested depths for 10~ 4 < v cut < 0.2 for both models of ionization 
energy loss. 

Results reported in this Section are evidence of accuracy in parameterizations for muon cross sections and sea level 
spectrum to be the principal source of uncertainties when simulating atmospheric muon flux at depths where neutrino 
telescopes are located. It contributes uncertainty from 3% (at the depth D — 1 km in pure water) to 15% (D — 
20 km) in resulting intensity of muon flux. Unfortunately, this level has at present to be considered as a limit for 
accuracy of muon propagation algorithms. Influence of model for ionization exceeds this limit only for mono-energetic 
muon beams with initial energies E < 10 TeV and only if level of observation is at very last stage of muon range 
where major fraction of initial muon energy has been lost. Actually, due to steep shape of atmospheric muon power 
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spectrum, an essential part of muons reaches detector location being just on the last part of path. Therefore effect 
remains noticeable also for real atmospheric muons but in this case uncertainty was found to be much less: 2-3%, 
which is in an excellent agreement with Refs. 12 23], while Ref. [^5| predicts much more significant difference (up 
to 17%). We suppose this disagreement may result from the fact that "small transfer grouping" technique used in 
Ref. 123] treats muon cross sections to be constant between two interactions in contrast with the MUM algorithm. In 
Ref. fp4f the same simplification was used but reported results were obtained by simulation with v cu t — 1CP 3 . With 
such small v cu t role of correct treatment for free path is not significant (see Sec. || and Fig. Q). Choice of value for v cu t 
is of even less importance and again, it is more critical if one investigates mono-energetic muon beam but for power 
spectrum alteration in v cu t within v cu t < 0.05 leads only to 1-2% differences in simulated muon flux intensities. Up 
to v cu t = 0.1 the error caused by rough account for fluctuations in energy losses remains less than one which comes 
from uncertainties with muon energy loss. This conclusion is in a good agreement with level of errors reported in 
Ref. [ p4[ . Differences between muon flux intensities simulated for different models of ionization and values of v cu t , as 
obtained in given work and in Refs. [f24|]25[], are presented in Fig. pi and Fig. ml. 
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FIG. 8. Intensity of vertical atmospheric muon flux I at different depth s D of pure water vs. v cu t as obtained by simulation 
with muons sampled according to sea level spectrum from Ref. |^] (Eq. ( ]4.l[ )). Closed circles: ionization is included in SEL; 
open circles: ionization is completely "continuous". Two horizontal solid lines on each plot show value for survival probability 
simulated with all muon cross sections multiplied by a factor k a = 1.01 (lower line) and k a = 0.99 (upper line) for v cut = 10~ 4 . 
Dashed lines on plots for D < 5 km correspond to intensity which was calculated for all energy losses treated as "continuous" . 
Dash-dotted lines show intensity of vertical muon flux simulated w ith ionization included in SEL, v cut = 10~ 4 and muons 
sampled according to the Gaisser sea level spectrum (Ref. |34j, Eq. (4.2)). Horizontal dotted lines correspond to v cu t = 10 -4 
and cross section for absorption of a real photon at photo-nuclear interaction parameterized according to Ref. ^] instead of 
parameterization proposed in Ref. P5(] which is basic in MUM. 
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FIG. 9. Dependencies for relation I2/I1 vs. water equivalent depth in standard rock as computed in this work (closed 
squares), in Ref. J24| (open squares) and in Ref. (closed circles). Ii is depth intensity for vertical atmospheric muon flux 
simulated with ionization included in SEL, I2 is one simulated with ent irely "continuous" ionization. Data for this work are 
obtained for sea level atmospheric muon spectrum from Ref. [3^1 (Eq. (4.1)) and v cut = 10 -3 ; data from Ref. represent 
result of simulation for sea level spectrum from Ref. j34j] (Eq. (4*2|)) and v cut = 10~ 3 ; data from Ref. (5^] were simulated with 
spectrum from Ref. fe*7| with "small transfer grouping" technique. 
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FIG. 10. Dependencies for relation I2/I1 vs. water equivalent depth in standard rock as computed in this work (closed 
squares) and in Ref. (open squares). Ii is depth intensity for vertical atmospheric muon flux simulated with entirely 
"continuous" ionization and v cut = 10 -3 , I2 is one simulated with the same treatment of ionization and v cut = 10 -2 . Data for 
this work are obtained for sea level spectrum from Ref. ||| (Eq. ([□])); data from Ref. H represent result of simulation for 
spectrum from Ref. @ (Eq. (Ejj)). 



So when simulating muon fluxes at large depths with an "ideal MC muon propagation algorithm" it is reasonable 
to use v cu t ~ 0.05-j-0.1 and entirely "continuous" model for ionization. Such setting of simulation parameters does 
not lead to the error which would be out of insuperable uncertainties with muon energy loss but allows to save the 
computation time essentially. Fig. 11 show dependence of computation time on v cu t and model for ionization, as was 
obtained with the MUM algorithm. Data for muon transportation codes PROPMU (Ref. §|]) and MUSIC (Ref. [||) 
are given on the figure, as well. We must emphasize that MUM in its presented version is ID algorithm, in contrast 
both to PROPMU and MUSIC. PROPMU treats only Coulomb multiple scattering while in MUSIC the angle of the 
muon acquired in all radiative processes is also simulated which takes an additional computation time. We evaluate 
the factor by which computation time with MUM would increase in case of extension up to 3D-algorithm as ~ 2. 
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FIG. 11. Averaged computation time T comp which is necessary for muon propagation in pure water vs. v C ut, as obtained 
with the MUM code. Thick lines correspond to muon with initial energy E a = 9 TeV transported down to D = 10 km. Thin 
lines are for E s = 1 TeV and D = 3 km. Solid lines show results for ionization included in SEL, dashed ones correspond to 
entirely "continuous" ionization. Circled asterisks on curves correspond to conservatively evaluated upper boundary for v cu t 
below which the MUM algorithm inner accuracy has been proved to be high enough. This limit is equal to v cut = 0.05 if 
ionization is included in SEL and to v C ut = 0.01 if ionization is entirely "continuous" (see Sec. Jni|). Circles and squares show 
values for T comp , as obtained with muon propagation codes PROPMU (version 2.01, 18/03/1993 with v cut = 10~ 2 which is 
unchangeable) and MUSIC (version for pure water with bremsstrahlung cross sections by Kehlner-Kokoulin-Petrukhin, 04/1999 
with v cut = 10 -2 which is unchangeable), correspondingly. Closed markers are for E s = 9 TeV and D — 10 km, open ones are 
for E s = 1 TeV and D = 3 km. 



Accounting for data on real accuracy of current version of the MUM code (see Sec. [II) which and data on compu- 
tation time presented in Fig. [ll] we conservatively consider v cu t = 0.05 and knock-on electron production included in 
SEL as an optimum setting for presented algorithm which allows to obtain the accurate results with relatively high 
speed. With such setting the proportion of computation time which is necessary to get the same statistics with MUM, 
PROPMU and MUSIC is approximately 1 : 10 : 600. Of course, for some methodical purposes it may be necessary 
to choose more fine v cut , e.g. if one wants to exclude an additional error when comparing results of simulations for 
different models of atmospheric muon sea level spectrum with each other or investigating survival probabilities which 
are much more sensitive to value of v cu t than simulated spectrum of atmospheric muons at large depths. 

We did not investigate specially the influence of simulation parameters on the results for the muon flux originated 
from neutrino but simple argumentation may be applied for this case. In contrast with atmospheric muons whose 
source is far away of underwater, under-ice or underground detector and whose flux may only decrease when passing 
from the sea level down to detector depth, the source for muons which are produced in vN interactions is uniformly 
distributed over water and/or rock both out- and inside the array. Intensity of the muon flux I® c which accompanies 
the neutrino flux in a medium is proportional to the muon range, and, consequently I® c cx [dE / dx)^Jj. al while 
simulated flux of atmospheric muons at large depths depends more sharply upon muon energy loss as was shown in 
this Section. Thus, one may conclude that the setting of parameters described above fits even better for propagation 
of muons originated from neutrino. 

It is impossible to foresee all particular cases and give some strict conformity between setting of parameters at 
muon MC propagation code and problem to be solved. But we tried to present in this Section the whole set of data 
which are necessary to chose the optimum set in each concrete case. 



V. SELECTED RESULTS AND COMPARISON WITH OTHER ALGORITHMS 



In this Section we present selected data on survival probabilities and atmospheric muon spectra deep underwater as 
simulated with MUM. To obtain atmospheric muon spectra we set v cut = 0.05. As was shown in Sec. Ill and Sec. IV 
it does not distort results comparing to simulation with smaller values of v cu t- To compute survival probabilities 
more delicate tuning was applied: v cu t = 10~ 3 . In both cases ionization was included in SEL. We compare our 
data with ones obtained with the PROPMU and MUSIC algorithms. Data simulated with PROPMU [version 2.01, 
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18/03/1993] (v cut = 10" 2 ) and PROPMU [version 2.1, 01/2000] (both with v cut = 10~ 3 and v cut = 10~ 2 ) are very 
close to each other, in all figures of this Section results from PROPMU [version 2.01, 18/03/1993] are presented. We 
used [version for pure water with bremsstrahlung cross sections by Kelner-Kokoulin-Petrukhin, 04/1999] with v cu t = 
10~ 3 for MUSIC. When comparing results on atmospheric muons at large depths obtained for pure and sea water the 
data are recalculated to each other using value p = 1.027 g cm - 3 as a sea water density (Ref. ]38| , |3"9|] ) . The difference 
between pure and sea water is negligible small for the muon propagation if one works in water equivalent units which 
was tested by us up to slant depth D = 10 km w.e. (see also Ref. |HJ). 

Fig. [12] shows survival probabilities vs. slant depth D in pure water as simulated for a set of initial muon beam 
energies from E s = 500 GeV to E s = 30 PeV. Survival probabilities obtained with MUM coincide within statistical 
errors with ones computed with MUSIC. PROPMU gives remarkably different values which are higher comparing to 
MUM and MUSIC output at muon energies E s < 30 TeV and become less at E s > 30 TeV. 




5 10 12.5 15 17.5 20 

D (km) 

FIG. 12. Survival probabilities vs. slant depth D in pure water as computed with MUM (solid lines), MUSIC (circles), and 
PROPMU (dashed lines). Figures near curves indicate initial energies of muon beams which were as follows: 500 GeV (1), 1 
TeV (2), 3 TeV (3), 10 TeV (4), 30 TeV (5), 100 TeV (6), 300 TeV (7), 1 PeV (8), 3 PeV (9), 10 PeV (10), 30 PeV (11). At 
simulations of data presented on the plot muons are treated as stopped as soon as their energy decreases down to 10 GeV. 



Fig. [13] gives more detailed data on survival probabilities for three particular cases. It presents muon spectra 
resulted from mono-energetic muon beams with initial energies E s = 1 TeV (Fig. |l3|(a)), E s — 9 TeV (Fig. |l3|(b)) 
and E s = 1 PeV (Fig. fl3|(c)) after propagation of distances 3 km, 10 km and 40 km in pure water, correspondingly. 
The distances were chosen so that survival probabilities would be much less than 1 in which case differences become 
more noticeable (see Sec. IV). A good agreement is observed between MUM and MUSIC data, while data obtained 
with PROPMU indicate essential differences which are of the same signs as in Fig. |l^. 

In Fig. ^ differential spectra for vertical atmospheric muons at different depths in pure water are present ed a s 
simulated with MUM, PROPMU and MUSIC. Muons at the surface were sampled according to spectrum Eq. (4.1). 
Also parameterizations for deep underwater muon spectra as proposed by A. Okada in Ref. |40|] and by S. Klimushin 
et al. in Ref. ]33| (from here on will call them "the Okada parameterization" and "the KBS parameterization", 
correspondingly) are shown. The KBS parameterization can adopt different models for sea level atmospheric muon 
spectrum. For data presented in Fig. [l4]we used spectrum Eq. (4.1) which is basic one for the KBS parameterization. 
MUM gives almost the same results as MUSIC which could be expected because survival probabilities for muons in 
pure water are the same when simulating with MUSIC and MUM, as was shown above. Simulation with PROPMU 
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produces the muon spectra which (a) are significantly higher (31%, 30%, 27% and 17% in terms of integral muon flux 
at the depths D = 1 km, 3 km, 6 km and 10 km, correspondingly) and (b) are expanded to the low energies. It is in 
good qualitative agreement with results on survival probabilities presented in Fig. [l^ and Fig. [l3| The coincidence 
between spectra simulated with MUM and curves for the basic KBS parameterization results from the fact that in 
both cases the same sea level atmospheric muon spectrum was adopted and, besides muon transport with the MUM 
algorithm was applied to obtain the KBS parameterization. We would like to mark that survival probabilities which 
KBS parameterization is based on were computed with v C ut — 10~ 3 . An excellent agreement with direct simulation in 
which v cut = 0.05 was set confirms the conclusion concerning insensitivity of results on simulated atmospheric muon 
spectra at large depths on value of v cut up to at least v cu t = 0.05 (see Sec. IV). The Okada parameterization is lower 
than KBS, MUM and MUSIC results (up to 18% in terms of integral muon flux at D = 1 km) at relatively shallow 
depths and becomes higher at D > 5 km because it is based on rather hard sea level atmospheric muon spectrum with 
index 7 = 2.57 (Ref. jyj) which leads to a deficit for low energy muons comparing to the basic KBS parameterization. 
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FIG. 13. Muon spectra resulting from mono-energetic muon beams with initial energies E 3 = 1 TeV (a), E s = 9 TeV (b) and 
E 3 = 1 EeV (c) after propagation down to depths D = 3 km, 10 km and 40 km of pure water, correspondingly, as simulated with 
MUM (histograms), PROPMU (circles), and MUSIC (triangles). Corresponding values for survival probabilities p (fraction 
of muons survived after propagation) are equal to: p(l TeV, 3 km) = 0.029 (MUM), 0.033 (MUSIC), 0.19 (PROPMU); p(9 
TeV, 10 km) = 0.030 (MUM), 0.031 (MUSIC), 0.048 (PROPMU); p(l PeV, 40 km) = 0.078 (MUM), 0.084 (MUSIC), 0.044 
(PROPMU). 



Fig. |l5| presents results on integral flux of vertical atmospheric mu ons at large depths in pure water as fa) simulated 
with MUM, PROPMU and MUSIC for sea level spectrum Ecl (O); (b) parameterized by KBS (Ref. Q) with sea 
level atmospheric muon spectra Eq. (4.1) (basic), from Ref. |Q (the Gaisser spectrum), from Ref. |l2| (the MACRO 
spectrum) and A. Okada (Ref. p0[) with sea level spectrum from Ref. plj ; (c) measured by S. Higashi et al. (Ref. p8|), 
V. M. Fedorov et al. (Ref. pjjj] ) and Yu. N. Vavilov et al. (Ref. Q). Note that "experimental" points on the plot does 
not represent the pure experimental data because authors had to recalculate obtained counting rates to the vertical 
direction using a model for the muon angular spectrum underwater. MUM and MUSIC results coincide with each 
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other within l-=-2%. Results from PROPMU algorithm exceed points from MUSIC and MUM by ~ 30% being higher 
than any of presented parameterizations. 

We also compared the data on muon propagation through the standard rock obtained with MUM and MUSIC. Mean 
energy for vertically down-going atmospheric muons sampled with sea level spectrum from Ref. |34f| was computed 
with MUM as£= 123±2 GeV, 256±4 GeV and 387±7 GeV at depths D = 1 km w.e., 3 km w.e. and 10 km w.e., 
respectively. The corresponding values simulated with the MUSIC code and reported in Ref. |24| are E = 125±1 
GeV, 259±3 GeV and 364±4 GeV. So the maximum difference observed at the depth D = 10 km w.e. is of 6%. 
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FIG. 14. Differential spectra of vertical atmospheric muons at four depths in the pure water as simulated 
PROPMU and MUSIC (in all cases muon energies at the sea level were sampled according to spectrum Eq. ([l.l|)) 
eterized according to KBS with sea level spectrum Eq. (fl.l|) (Ref. |B3[) and A. Okada (Ref. |4Q|). 



with MUM, 
and param- 



Thus, results on survival probabilities and atmospheric muon spectra at large depths as simulated with MUM are 
practically in a coincidence with ones obtained with MUSIC and are not also in contradiction with any experimental 
and theoretical results presented in this section. The PROPMU algorithm shows noticeable differences with MUM 
which are in good qualitative agreement to each other: higher survival probabilities lead to higher muon fluxes deep 
underwater. It is difficult to clarify the source of observed discrepancies without detailed comparison for all steps of 
the algorithms but we believe that they can not be explained only by a difference in models for muon energy loss as it 
is used in MUM (see Appendix |X|) and PROPMU (Refs. which does not exceed 2% at E < 10 TeV (in terms 

of stopping power) being besides of both signs. 

More data obtained with the MUM algorithm can be found in Ref. |B3| . 
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FIG. 15. Results for integral flux of vertical atmospheric muons vs. depth in pure water as (1) simulated with MUM, 
PROPMU and MUSIC with sea level spectrum Eq. ([uj]); (2) parameterized by KBS (Ref. Q) with sea level atmospheric 
muon spectra Eq. ( [O] ) (basic), from Ref. El (the Gaisser spectrum), from Ref. E^] (the MACRO spectrum) and A. Okada 
(Ref. Q) with sea level spectrum from Ref^fuJ ; (3) measured by S. Higashi et al. (Ref. ||), V. M. Fedorov et al. (Ref. Q) 
and Yu. N. Vavilov et al. (Ref. El). 



VI. CONCLUSIONS 



We have presented the muon propagation Monte Carlo FORTRAN code MUM (MUons+Medium) and have given 
selected results obtained with the code for muon spectra at large depths and survival probabilities in comparison with 
results obtained with other muon transportation algorithms. It was shown that for majority of applications it is quite 
enough to account only for fluctuations in the radiative energy loss with fractions of energy lost as large as v > v cu t — 
0.05-^0.1 while ionization energy loss may be entirely accounted by stopping power formula, as well as radiative energy 
loss with fractions of energy lost v < v cut = 0.05^-0.1. This gives an essential advantage in terms of computation 
time comparing to commonly used v cu t — 1(U 3 -^10 -2 without lost of accuracy when simulating both propagation of 
atmospheric muons and muons which are born in uN interactions. However in practice it makes particular demands 
to accuracy of MC algorithm. Some customary simplifications (e.g. Eq (2.4)) which work perfectly when v cu t = 
1CU 3 -^10 -2 become sources of significant errors when v cut increases. The sign and value of these errors depend also 
on whether fluctuations in ionization are accounted or not. So for presented version of the MUM algorithm the 
optimum set of simulation parameters was conservatively evaluated by us (accounting results on inner accuracy test 
and dependence of computation time on v cu t) to be v cu t = 0.05 and knock-electron production included in SEL. 

Our point of view on advantages of MUM is as follows. It is flexible enough and provides with eventuality to tune 
parameters of simulation to an optimum for each concrete case to get desirable equilibrium between computation 
time and accuracy. Medium composition and parameterizations for muon cross sections are easily changeable. Inner 
accuracy of the code was conservatively evaluated to be 2xl0~ 3 or better. Besides, MUM provides with the special 
routine which allows to test inner accuracy for each given set of simulation parameters and take it into account when 
evaluating significance of results. The main disadvantage is that MUM in its reported version does not still treat 
the three dimensions like, e.g. PROPMU and MUSIC codes. So it can not be used to obtain lateral and angular 
deviations of muons at propagation through matter. Also other important features are still missed in MUM - for 
instance, treatment of composed medium as it is possible, e.g. in the latest version of the MUSIC code (Ref. iflBI). 
But we consider the current version of MUM as a basis for further development and plan to complement it, step by 
step, by all necessary features. 

The MUM code is available on request to sokalski@pcbailO.inr.ruhep.ru. 
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APPENDIX A: PARAMETERIZATIONS FOR MUON CROSS SECTIONS USED IN THE MUM 

ALGORITHM 

We use the following designations in this section: a — 7.297353x 10 -3 - fine structure constant; r e = 2.817941 x 10 -13 
cm - classical radii of electron; m M = 0.1056593 GeV and m e = 0.5110034 MeV - muon and electron masses, corre- 
spondingly; Na — 6.02 2 x 10 23 - the Avogadro number; Z and A - electric charge and atomic weight, correspondingly; 
e = 2.718282; ir = 3.141593. Other notations are explained in comments to formulas when necessary. 



1. Bremsstrahlung 

We use differential cross section for bremsstrahlung as parameterized by Yu. M. Andreev, L. B. Bezrukov and E. 
V. Bugaev in Ref. E9] as a basic parameterization: 
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Integration limits for bremsstrahlung in Eqs. (2.1), ( |2.3|) , (2.6) and (3.3) arc 
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»L = 0, v b max = 1 - -^{m^/E)Z 



1/3 



Note that this parameterization does not account for contribution from e-diagrams for bremsstrahlung when 7- 
quantum is emitted by atomic electrons which are knocked on by recoil (Ref. |j47| ) . Corresponding corrections are 
done in parameterizations for knock-on electron production (see Appendix A 4) according to Ref. E§]. Optionally 
differential cross section for muon bremsstrahlung can be also treated in MUM according to parameterization given 
by S. R. Kelner, R. P. Kokoulin, and A. A. Petrukhin (Ref. MjM). 



2. Photo- nuclear interaction 

We use parameterization for the photo-nuclear interaction of muon proposed by L. B. Bezrukov and E. V. Bugaev 
(Ref. H): 
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z = 0.00282A 1/3 cr 7A r, t = Y if — , ml = 0.54 GeV 2 , m\ = 1.80 GeV 2 . 

Total cross section for absorption of a real photon of energy v = s/2m^ = vE by a nucleon, c 7 Ar, can be calculated 
in MUM optionally according to either parameterization from Ref ]35| (basic) : 

<7yN = [114.3+ 1.6471n 2 (0.0213i/)]^6, 

or by the ZEUS parameterization (Ref. |}(| ) : 

a lN = (63.5 s - 097 + 145 s" ' 5 ) nb, 

where s and v are expressed in GeV 2 and GeV, correspondingly. 

Parameteriz atio n (Ref. p5|) is vali d for v > 1 GeV so we use values !)™ in = 0.8/-E(GeV) and = 1 as integration 
limits in Eqs. (2A_), (t?-3j), ( P^q) and (3J3). Note that results of integration were tested to be almost insensitive to the 
lower limit in a wide range 0.2/S(GeV) < w" m < l.5/E(GeV). 



3. Direct electron-positron pair production 



Cross section for direct e + e -pair production is used in MUM as parameterized by R. P. Kokoulin and A. A. 
Petrukhin (Refs. MM)- 
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Here p = (e + — e )/(e + + e ) is the asymmetry coefficient of the energy distribution of e + e -pair, e + and e arc 
positron and electron energies, correspondingly. Limits for integration over p are determined by: 



0<|p|< 1 



E 2 (1 - v) 
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4m e 
Eu 



i? is a parameter determined by the value of radiation logarithm (R = 183 for Thomas-Fermi model and slightly 
depends upon Z for Hartrey-Fock model). Its values are taken from Ref. |5C|j , where R has been calculated for 
different atoms according to Hartrey-Fock model. (,(Z) ~ 1 takes into account the pair production in coll isio ns wit h 
electrons. Values from Q(Z) are computed according to Refs. [^8|,^l). Integration limits for j = p in Eqs. (2.1), (|2+|), 
(2.6) and (3.3) are 
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4. Knock-on electron production 



Cross section for knock-on electron production is parameterized in the MUM algorithm as follows: 



da 
dv 



^-(E,v) = 2nr 2 Z^l---. 



1 1 E 



\ ) + A, -IE. 



2m P E 



m 2 + 2m e E 



A ej (E,v) represents the correction which takes into account e-diagrams for bremsstrahlung (Refs. [ |47|j4q ] ) resulting 
in additional recoil electrons: 
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value of d^j. is used also as upper integration limit in Eqs. (2.3) and (2.6) for j = e. 
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5. Ionization 



Following Refs. we treat in the MUM code e-diagrams for bremsstrahlung as a part of ionization process. 

Therefore we have to use a bit modified formula for ionization: 
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Here K = 0.1535 MeVg 1 cm 2 , p is the muon momentum, (3 = p/E is the muon velocity, p is the material density, / 
is the mean ionization potential, 



E„ 



(2m e p 2 )/(m 2 + m 2 e + 2m e E) 



is the maximum energy transferable to an electron, S is the density-effect correction which is treated according to 
Ref. [El: 



5 = 6(X - X ) [4.6052X + a6(X 1 - X)(X 1 - X) 



C] 



where is the step function (9(x) = at x < and 6(x) = 1 at x > 0), X = log 10 (p/m M ). The values Xq, X\, a, 
m and C depend on the material and can be found in Refs. Jh],|5^] along with values for J, p and Z/A. The first 
term represents Bethe-Bloch formula with corrections for density eff ect, the second one accounts bremsstrahlung e- 
diagrams. Expressions for A e7 (E,v) and i^ a;c are given in Appendix A 4, meaning of values A e ff and ki is explained 
in Sec. || 



APPENDIX B: FREE PATH BETWEEN TWO MUON INTERACTIONS 



For the proof of the set of equation Eqs. ( |2.2| ) it is convenient to introduce the kinetic equation for a propagation 
of a mono-energetic muon beam through a medium. With the notations used in textbooks this equation has the 
following view: 



dn{E,t)/dt - d[/3(E)n(E,t)}/dE + n(E,t)/\{E) = 

n(E,0) = n 5{E-E a ) 



(Bl) 



Here, n(E, t) is the number of muons with energy E after propagation of distance t, (3(E) is the "continuous" energy 
loss per unit path, X(E) is the muon mean free path before interaction of SEL type. The solution of Eq. (Bl) is: 



n{E, t) = n 5{E - e(E ,t)) exp 
where e(Eo,t) is found from the equation 
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can be rewritten as: 
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(B4) 



where rj(E,t) is the pro babi lity fo r a single muon to pa ss t he path t without interaction of SEL type and then one 
can easily see that Eqs. (B3) an d (p34f) lead to the Eqs. (2.2) after the following substitutions which are necessary for 
a return to the notations of Sec. O. 

E E -> E 2 , \{E) -» L(E), (3(E) -» [dE(E)/dx] CEL , t -» L. 
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